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PREFACE 


This  report  was  prepared  by  Dr.  Robert  B.  Herrmann  of  St.  Louis 
University,  St.  Louis,  Missouri,  as  part  of  ongoing  work  at  the  U.  S. 
Army  Engineer  Waterways  Experiment  Station  (WES)  in  Civil  Works  Invest! 
gations,  "Methodologies  for  Selecting  Design  Earthquakes,"  sponsored  by 
the  Office,  Chief  of  Engineers,  U.  S.  Army. 

Preparation  of  the  report  was  under  the  direction  of  Dr.  Ellis  L. 
Krinitzsky,  Chief,  Engineering  Geology  Research  Facility,  WES.  General 
direction  was  by  Mr.  James  P.  Sale,  Chief,  Soils  and  Pavements  Labora- 
tory, and  Mr.  Don  C.  Banks,  Chief,  Engineering  Geology  and  Rock  Me- 
chanics Division. 

COL  John  L.  Cannon,  CE,  was  Commander  and  Director  of  WES; 

Mr.  F.  R.  Brown  was  Technical  Director. 


2 


1a 


iBiw'»aat 


c 

[ 

( 


1 1 
*• 

f 


CONTENTS 


PREFACE  

PART  I:  INTRODUCTION  

PART  II:  THEORY 

Statement  of  tne  Problem 

PART  III:  NUMERICAL  TECHNIQUES  AND  TESTS  . . 

Contour  Integration  .....  

Numerical  Integration  Techniques.  . . . 
Independent  Solution  for  Simplified  Case 
Numerical  Tests  

PART  IV:  SUMMARY  AND  CONCLUSIONS  

REFERENCES 

APPENDIX  A:  COMPUTER  PROGRAM  WESHASK  .... 

Function 

Input  

Output 

File  Codes 

Sample  Input 

Program  Listing  

APPENDIX  B:  COMPUTER  PROGRAM  WESREFL  .... 

Function 

Input  

Output 

Sample  Input  

Program  Listing  


Page 


2 

4 

5 

5 

12 

12 

16 

19 

24 

4b 

50 

A1 

A1 

A1 

A1 

A1 

A2 

A2 

B1 


*. 

1 


B1 

B1 

B1 

B1 

B1 


k 


NOTE:  Appendixes  A and  B have  been  deleted  to  comply  with 
ER  18-1-6. 


3 


EARTHQUAKE  GENERATED  SH  WAVES  IN  THE  NEAR  FIELD 
AND  NEAR-REGIONAL  FIELD 


PART  I:  INTRODUCTION 

1.  The  description  of  earthquake  generated  ground  motion  is  a 
difficult  problem  due  to  the  complex  nature  of  the  earthquake  process 
itself  as  well  as  the  process  of  transmission  of  seismic  energy  through 
the  real  earth.  Because  of  this  complexity,  it  is  difficult  to  fit  ob- 
served earthquake  strong  motion  data  to  simple  prediction  models  with 
any  reliability.  Many  attempts  have  been  made  to  fit  strong  motion  data 
as  a function  of  parameters  such  as  earthquake  magnitude,  focal  depth 
and  epicentral  distance.  Such  models  usually  succeed  in  predicting  the 
mean  of  the  observed  data  set,  but  because  of  their  strictly  empirical 
nature  offer  no  insight  to  the  cause  of  the  scatter. 

2.  Some  of  the  error  in  the  fit  of  these  predictive  models  can  be 
accounted  for  by  making  the  model  more  complex.  However,  this  cannot  be 
done  blindly  since  the  addition  of  the  extra  parameters  required  to  ex- 
plain the  scatter  of  simpler  models  may  not  have  a theoretical  justifi- 
cation. In  any  case,  a model  based  upon  a strictly  empirical  fit  may 
not  be  general  enough  to  admit  Editions  to  the  data  base.  A better  ap- 
proach to  the  modeling  of  obs  data  is  to  perform  a theoretical 
study  of  the  generation  of  ".qu^ke  ground  motion  so  that  an  insight 
on  the  choice  of  the  empirical  model  used  can  be  obtained. 

3.  The  purpose  of  this  report  is  to  present  work  performed  on 
modeling  the  SH  component  of  ground  motion  for  elementary  earthquake 
point  sources  for  ground  motions  with  frequencies  less  than  1.0  Hz  at 
observation  sites  5-500  km  from  the  source.  While  the  extension  of  this 
report  to  include  higher  frequencies  is  important  for  the  understanding 
of  theoretical  SH  wave  propagation  in  these  distances  ranges,  the  spa- 
tial inhomogeneities  of  the  real  earth  are  such  that  an  analytical  solu- 
tion is  not  obtainable  at  frequencies  much  greater  than  1 Hz  for  a real 
earth.  However  the  insights  of  the  modeling  performed  here  will  be  of 
great  value  in  understanding  the  gross  properties  of  the  high  frequency 
content  of  earthquake  generated  motion. 
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PART  II:  THEORY 


4.  A recent  paper  by  Helmberger  and  Malone^  demonstrates 

how  well  earthquake  ground  motion  can  be  modeled  by  taking  into 

account  the  effect  of  crustal  layering  on  the  seismic  signal.  They 

2 3 

used  a generalized  ray  technique  * to  construct  their  solution. 

This  technique  involves  the  computation  of  all  possible  contributions 
to  the  seismic  signal,  whether  due  to  multiple  reflection  or  refrac- 
tion. The  advantage  of  this  technique  is  that  it  is  valid  for  high 
frequencies.  A disadvantage,  though,  is  that  for  a reasonable  model 
of  the  earth  the  number  of  individual  ray  contributions  to  the  seismic 
signal  becomes  so  large  that  the  computer  algorithm  required  is  also 
very  complex.  The  generalized  ray  technique  is  thus  limited  to 
either  simple  earth  models  or  to  studies  of  just  part  of  the  total 
ground  motion  time  history. 

4 5 

5.  Another  approach  is  that  developed  by  Haskell  and  Hudson. 
This  method  does  not  consider  the  contribution  of  the  individual 
reflected  and  refracted  arrivals,  but  rather  yields  the  complete 
solution  for  a source  in  a layered  halfspace.  The  drawback  of  the 
wave  theory  approach  is  the  difficulty  of  finding  solutions  at  high 
frequencies.  However,  the  earth  model  can  be  as  complicated  as  re- 
quired. It  is  the  wave  theory  method  that  is  the  subject  of  this 

report.  The  basic  theory  used  is  a modification  of  that  developed  by 
4 5 

Haskell  and  Hudson.  The  modifications  involved  recasting  their 
solutions  into  forms  amenable  to  the  numerical  computation  of  not  only 
the  surface  wave  but  also  the  refracted  and  reflected  body  wave  con- 
tributions. 


Statement  of  the  Problem 

6.  The  problem  considered  is  that  of  SH  wave  generation  by  an 

a.  h 

elementary  point  source  in  the  m layer  of  an  N layered  structure 
with  an  upper  free  surface.  Each  layer  is  homogeneous  and  isotropic 


with  compressional  wave  velocity  , shear  wave  velocity  6^,  and 
density  (k*  1,  N).  The  Ntfl  layer  extends  dov.’nwards  to  infinity. 

A cylindrical  coordinate  system  (r,<|>»z)  is  used  with  origin  on  the 
free  surface  directly  above  the  source,  with  the  z-axis  taken  positive 
downwards.  The  layer  interfaces  are  the  planes  z = zk  (k  = 1 ,2,...,N-1 ). 
The  source  is  situated  on  the  plane  z » zm_j  + hm.  For  purposes  of 
derivation,  the  source  is  restricted  to  lie  in  one  of  the  layers  above 


th 


the  halfspace.  (T  <=  zm  - zm Is  the  thickness  of  the  m layer. 

7.  Hudsonb  gives  the  following  expression  for  the  Fourier  trans- 
form of  the  displacements  at  the  free  surface  z ■ 0: 


u,(r,4>,  | cos  n<p  +g”  sin  np)  '4' 

■ *0  0 

- -J-  (*/*  cos  np  +*/*  sin  np)  jm(kr)/FL  J 

m m 

u+(r,4>,  a,  o>)  = ^ j dk  j~ (£,** cos np -g* sin np)  Jm(kr)/F  (1 ) 

- (g“ cos np-gf* sin n<f>) J 

m m 

u,(r,  p,  o,  c a)  *=  cos  np  +g,u  sin  nip)  Jm(kr)/F  J . 

In  this  expression  u^,  the  radial  displacement,  is  positive  in  a direc- 
tion away  from  the  source;  u^,  the  tangential  displacement,  is  positive 
in  a direction  of  increasing  tTz,  the  vertical  displacement,  is  posi- 
tive downward.  The  explicit  expressions  for  the  various  terms  in  these 
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equations  are  given  by  Haskell  and  Hudson. 

8.  Care  must  be  taken  In  evaluating  the  expressions  in  Equation  1 
because  of  the  presence  of  non-causal  terms. ^ The  non-causality  of 
certain  parts  of  the  solution  arises  because  cylindrical  potentials 
were  used  to  form  the  solution,  when  a cartesian  coordinate  system 
v/ould  be  more  appropriate.  Because  of  this,  the  nature  of  the  solution 
of  Equation  1 is  such  that  very  close  to  the  source  the  tangential  and 
radial  displacements  contain  P,  SV,  and  SH  wave  components. 
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9.  The  solution  of  Equation  1 is  simplified  if  one  is  interested 
primarily  in  the  propagating  wave  contribution.  Using  the  relation 

|p  Jn(kr)  = k Jn-1(kr)  - J Jn(kr)  (2) 

and  ignoring  terms  decreasing  faster  than  r~\  the  solution  for 
the  Fourier  transform  of  the  tangential  displacements  becomes 


u,(r,4>,0,<o)  = l / -(g"s  cos  n$ 
♦ n=0  0 ^ 


gjc  sin  n*)Jn.1(kr)F[1k  dk  (3) 


where 


gnC  = 
y* 

(L21 

-L„) 

+ <42  - 

L ) s'1C 

L12'  s2 

gnS  = 

(L21 

-L„) 

*?S  * (l22  - 

L 1 sns 
l12'  s2 

(4) 

FL  “ 

J11 

' J21 

• 

10.  The  and  E^  are  elements  of  the  L and  E matrices  which 
are  defined  by  the  matrix  products 

J = £f/  AN-1  ^dN-l  ^ — A1  ^dl  ^ 

(5) 


l*  en] 


The  various  matrices  in  Equation  5 are  defined  as 


pnV 


0 


0 


(6) 


7 


(7) 


A(z) 


and 


0(2)  = 


V‘>B\ 


p&  Ve 

e.  1 

V* 

-VpvS 

2 

ft  v S 
ee 

2 

'e  ce 

Sg(z)  * sinh 

VgZ,  and 

(8) 


The  elements  of  the  matrices  are  to  be  evaluated  using  the  layer 
parameters  of  the  particular  indicated  by  the  subscripts  of  Equation  5. 

11.  Tne  source  coefficients  SjC,s  are  adapted  from  Haskell.^ 

It  is  generally  accepted  that  a system  of  point  forces  oriented  to 
form  a double-couple  without  moment,  or  equivalently  two  perpendicular 
dipoles  of  opposite  sign,  is  an  adequate  representation  of  the  dis- 
location model  of  an  earthquake  source.  Movement  on  a fault  plane 
can  be  described  by  either  of  the  two  force  systems  mentioned  above 
or  by  defining  the  direction  of  slip  on  a fault  plane  of  given  dip 
and  strike.  The  interrelationship  of  these  three  ways  of  describing 

7 

the  earthquake  source  is  given  by  Herrmann.  The  use  of  the  perpen- 
dicular dipole  representation  leads  to  the  least  confusion  in  applica- 
tion. Let  the  orientation  of  the  tension  and  pressure  axes  in  car- 
tesian coordinates  be  given  by  the  vectors  T - (fpfg.fg)  and 
n * (n^^.n^),  respectively.  The  source  coefficients  s!?c  and  s^ 
are  all  zero  except  for  the  terms 


S!C  = - 2<flf3  - "l"3>/4*4 
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- 2<f2f3  • n2n3)/4*Bm 


sfC  = k(f2  “ fl  ' n2  + nl)/4,T4  (9) 

s2S  * - 2k(f ^ f2  - n1n2)/4irB^ 

where  k is  the  wavenumber  and  is  the  shear  wave  velocity  in  the 

source  layer.  The  expressions  for  , D,  and  svc’s  differ  from  those 
4 5 w 1 

given  by  Haskell  and  Hudson  in  that  they  have  been  modified  to  elim- 
inate apparent  numerical  singularities.  The  final  solution  given  by 
Equation  3 is  the  same  as  that  given  by  these  authors. 

12.  For  an  earthquake  source  described  by  the  source  coefficients 
of  Equation  9,  the  desired  solution  of  Equation  3 can  be  rewritten  as 
tne  sum  of  two  terms: 


= F1(f,n,(j))G1  + F2(f  ,n,<t>)G2 


where 


/ -(91/FL)J0(kr)k  dk 
0 

l -(92/FL)J1(kr)k  dk 


L21  ' L11 


g2  ~ k^L22  ” L12^ 


°1 1 ' °21 


2 

F-j (T,n ,4> ) = - 2[(f2f3  - n2n3)cos^  - (f-j f3  - n^njjsin^  ]/4nBm 

F2(f,n,<j>)  = - C2(f  1 f2  - n^Jcos  2*  + (f2-f^-n2+n^)sin  2<j1]/4ne*. 
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13.  To  understand  the  source  terms,  three  simple  models  of 

an  earthquake  are  considered.  Let  the  x and  y axes  of  the  cartesian 
coordinate  system  point  north  and  east,  respectively.  The  angle  t 
is  then  the  azimuth  measured  clockwise  from  north. 

i.  Right  lateral  vertical  strike-slip  fault  striking  north 

f = (-.707, .707,0}  n=  (.707, .707,0) 

= 0 

2 

F2(f,n,<t»)  = 2 cos  2$/4*&m 

b.  Thrust  faulting  on  a fault  dipping  45°  to  the  east  or 
west  and  striking  north  (45°  dip  slip) 

f-  (0,0,1)  n=  (0,1,0) 

F1  (?’,n,4>)  = 0 

2 

F2(f,n,$)  = sin  2<fr/4*3m 

c.  Vertical  dip  slip  faulting  on  a fault  striking  north 
with  the  east  side  downthrown 

f = (0,. 707, .707)  n"  = (0,-. 707, .707) 

— _ 2 
^(f.n,*)  * - 2 cos$/4TTBm 

F2(f,n,$)  = 0 . 

14.  An  examination  of  these  terms  indicates  that  for  vertical 
strike  slip  faulting  or  for  45°  dip  slip  faulting,  the  solution  for 
SH  wave  motion  is  proportional  to  the  function  G2  only.  Likewise  for 
pure  dip  slip  on  a vertical  fault,  the  SH  wave  motion  is  a function  of 
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G-j  only.  It  Is  also  Interesting  to  note  that  a vertical  strike  slip 
fault  is  twice  as  efficient  a generator  of  SH  waves  as  a 45°  dip  slip 
fault.  The  numerical  techniques  and  examples  to  follow  deal  with  the 
functions  Gj  and  Gg  only,  and  must  be  multiplied  by  F^(f,n,«|>)  and 
f°r  proper  scaling  before  the  synthetic  seismograms  are 
compared  to  real  data. 

15.  After  Equation  10  is  evaluated,  it  must  be  multiplied  by 
the  Fourier  transform  S(w)  of  s(t),  the  source  time  function.  For  a 
dislocation  source,  s(t)  represents  the  time  history  of  the  faulting 
process  for  which  s(t)  = 0 for  t < 0 and  s(t)  - fiQ  for  t » 0.  The 

seismic  moment  MQ  is  defined  by  the  relation  MQ  * uuA,  where  jj  is  the 

rigidity  of  the  medium  in  which  the  faulting  occurs,  U is  the  average 
fault  displacement  and  A is  the  fault  area.  In  normal  use  the  CGS 

units  of  Mq  are  dyne-cm.  The  SH  wave  ground  motion  as  a function  of 

time  is  obtained  by  taking  the  inverse  Fourier  transform  of  the 
product  of  S(u  ) and  Equation  3: 

l * _ 

u (r,<f>,0,t)  = (2n)  / $(w)u  (r,<j>,0,a>)exp(ia>t)  du  . (12) 

” -00  ™ 

5 

Hudson  discusses  the  mathematic^  of  superposition  of  solutions  such 
as  Equation  12  to  represent  extended  sources  and  complex  rupturing 
processes.  However,  these  extensions  are  beyond  the  scope  of  this 
report. 
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PART  III:  NUMERICAL  TECHNIQUES  ANO  TESTS 


16’  As  seen  In  the  previous  section,  the  solution  of  SH  wave 
motion  due  to  an  arbitrary  point  dislocation  source  model  of  an 
earthquake  reduced  to  the  evaluation  of  two  Integrals, 


Gt  ■ /'(91/FL)J0(kr)k  dk 

0 

G2  * /“(92/FL)Ji(kr)k  dk 


(13) 


where  g-j,  g2,  and  FL  are  defined  In  Equation  11.  The  objective  of  this 
section  Is  to  recast  these  Integrals  Into  a form  suitable  for  numerical 
Integration.  As  an  Independent  test  of  the  numerical  Integration 
technique,  the  problem  is  solved  by  a different  technique  for  a simple 
earth  model,  and  the  two  solutions  are  compared. 

Contour  Integration 


17.  A contour  Integration  method  for  evaluation  of  integrals 
of  the  form  of  Equation  13  has  been  described  by  many  authors.  The 
development  given  here  parallels  one  presented  by  Ewing,  Jardetzky 
and  Press. G 

18.  The  two  integrals  in  Equation  13  are  of  the  general  form 

I = / f(k,0))Jn(kr)k  dk  . (14) 

0 

The  function  f(k,u)  has  a finite  number  of  poles  on  the  real  axis 

for  k < k < k and  a branch  point  at  k = k . From  Equation  13, 
bN  ~ “ H eN 

f(k,dj)  * 9]/fl  For  n * 0 and  f(k,&>)  ■ g2/F^  for  n ■ 1.  A study  of  the 
functional  forms  of  g^,  g2,  and  shows  that  f(k,w)  is  even  In  k for 
n even  and  odd  in  k for  n odd. 

19.  The  Bessel  function  of  the  first  kind  can  be  expressed  in 
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terms  of  the  Hankel  functions  of  the  first  and  second  kind  by  the  relation 
Jn(kr)  - ]■  [ H^^kr)  + H<2)(kr)  ] . (15) 

Using  Equations  14  and  15,  one  obtains 


1*2-/  ffk.ujH^^krJk  dk 
o 

+ \ 7 f(k,o>)H<2>(kr)k  dk 


(16) 


20.  The  integrals  of  Equation  16  are  improper  since  the 
integrands  become  infinite  at  the  zeros  of  FL.  To  evaluate  these 
integrals,  contour  integration  is  performed  in  the  complex  k plane. 
The  contours  for  evaluating  and  Ig  are  shown  in  Figures  la  and  lb, 
respectively. 


Fig  1.  Contours  in  complex  k plane  for  evaluating 
the  integrals  1^  and  ^ 


In  these  figures,  A is  the  branch  point,  P is  a pole  on  the  real  axis, 
and  the  heavy  line  along  part  of  the  real  axis  and  along  the  negative 
imaginary  axis  is  the  branch  cut.  These  contours  are  chosen  since 
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H^(kr)  is  analytic  in  the  first  quadrant  while  H^(kr)  is  analytic 
in  the  fourth  quadrant  of  the  complex  k plane.  In  the  first  quadrant. 
Ini  v.  >0  and  Re  va  > 0,  while  in  the  fourth  quadrant  Im  \>a  <0  and 

sN  - sN  - eN  - 

Re  vQ  > 0. 

eN  - 

21.  Using  Cauchy's  integral  theorem.  Equation  16  becomes 

I = \ I f(c.u)H^(cr)c  d 5 + i / fU,w)H*;2*(crk  dc 

c 0 n 2 E n 

+ j / f(t,o>)H*2*(cr)c  dc  - *1[Res[f(c,u)H^2^(cr)c  (17) 

2 OAE  n n 


+ \ / f(5»u)H^(cr)c  dc  + J-  / f(c,u)H^2^(cr)c  dc  . 

1 r2 


As  the  arcs  and  r2  extend  outward  to  Infinity,  the  contribution 
of  the  last  two  terms  in  Equation  17  to  the  total  solution  goes  to 
zero.  Using  an  expression  for  the  analytic  continuation  between  the 
Hankel  functions 

H^1}(cr)  * - exp(-1n*)Hj2)(-cr), 

Equation  (17)  becomes 


I 


2"  / [f_U,<*>)  - exp(-inw)f+(-c,w)]H^(cr)c  dc 

o 

kg 

+ i / [f.(k,w)  - f (k,u)]H^2)(kr)k  dk 
c 0 " n 


(18) 


- *1  I Res[f(c,u)H^2^(cr)c]  , 


where  the  + or  - subscripts  Indicate  that  Im  v.  > 0 or  Im  v.  < 0, 

eN  eN 

respectively,  is  to  be  used  to  evaluate  the  expressions  for  f(c,u>). 
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22.  By  a simple  change  In  variable  Equation  (18)  can  be 
rewritten  as 

1 * " 2"  o - exp(-1nu)f+(lT  ,u)]H^(-1t)t  dT 

ke 

+ l / [f+(k,<o)  - f (k ,o>)]H^(kr)k  dk  (19) 

- *i  l Res[f(c,u)H^(cr)c]  . 

23.  The  final  steps  of  the  derivation  make  use  of  specific 
properties  of  the  functions  g^,  gg,  and  F^.  Letting  f*(c,w)  * g-j/F^ 

and  f^(j;,w)  * g2/FL,  **  can  be  s^10wn 
f|(x,u)  - f}(-X,a)) 

fj(x,u)  * - fj(-x,w) 

f^(-ix,u)  * f|( 1x7u ) (20) 


f^(-1x,<o)  = f*(1x,u>) 


fl’2(x,<o)  - fl»2(x,u.)  , 

where  x is  a real  quantity  and  the  bar  above  the  function  Indicates 
that  the  complex  conjugate  should  be  taken. 
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1 


24.  The  modified  expressions  for  Gj  and  G2  of  Equation  13 
finally  become 

G,  = i / Im[f}(k,ui)]H^2)(kr)k  dk 


- *i  l ResCf^c.ujH^^cr)?] 


" f o Im^+^T»“)3K0(Tr)T  dT 


(21a) 


Go  = 


bn 

i / Im[f2(k,u)]H!2)(kr)k  dk 

0 1 


- *1  l Res[f2(c,w)H{2)(cr)c] 


+ - / Re[f+(iT,u)]K1(ir)T  dx 


where  the  additional  relations 

2 


and 


T7  Vz> 


Hj2)(-iz)  = | K,(z) 


(21b) 


have  been  used.  Kn(z)  is  a modified  Bessel  function  of  order  n. 


Numerical  Integration  Techniques 

25.  A computer  program,  WESHASK,  has  been  written  to  perform  the 
computations  required  to  evaluate  Equation  21.  Appendix  A contains  a 
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listing  and  a description  of  the  use  of  this  program.  It  is  approprl 
ate  to  discuss  the  numerical  integration  techniques  used  as  well  as 
their  limitations.  Each  of  the  three  terms  making  up  and  G2  in 
Equation  21  is  discussed  separately. 

26.  The  first  term  in  both  G^  and  Gg  Is  an  integral  along  the 
real  axis  of  the  form 


f(k,u>)H^(kr)k  dk 


(22) 


To  evaluate  this  integral,  a change  of  variable  Is  introduced.  Using 
relation  k * k.  sin  y , Equation  22  becomes 

ir/2  /n\  o 

/ f(k  siny,o))H^'(k  r siny)k£  siny  cosy  dy  (23) 

0 eN  n tiN  eN 


The  advantage  of  Equation  23  is  that  the  Integrand  has  a zero  weight 

at  k * kD  , where  the  branch  point  may  introduce  a singularity  at 
6N 

some  frequencies,  as  well  as  at  k * 0,  where  the  Hankel  function  be- 
comes undefined. 

27.  Equation  23  is  evaluated  by  using  a trapezoidal  integration 

rule: 


j f(kj »“)H^(kir)k^siny.  cosYi  Ay  , 

where 

Ay  * ir/2M  , y.  .=  iAy  , andk.  =k  siny.  . 

1 1 8N  1 


(24) 


28.  The  choice  on  the  summation  index  M depends  upon  both 
the  maximum  angular  frequency,  n>,  and  epicentral  distance,  r,  consi- 
dered. The  Hankel  function  H^(z)  Is  an  oscillatory  function.  To 
avoid  numerical  Integration  problems,  the  Hankel  function  must  be 
sampled  at  least  twice  per  oscillation.  For  a given  frequency  f 
(w*  2irf)  and  distance  r,  the  value  of  M should  be  large  enough  to 
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MBS* 


"TrniMrir.rar.ai 


to  approximate  the  inequality 


M > nfr/eN  . 

Good  numerical  results  have  been  obtained  using  M * 100  for  frequencies 

less  than  O.b  Hz  for  distances  less  than  500  km  with  6^  = 4.67  km/sec. 

Even  though  tnis  value  of  M does  not  satisfy  the  inequality,  it  is 

reasonable  to  expect  that  the  guidance  given  by  the  above  inequality 

yields  the  proper  order  of  magnitude  for  M. 

2S.  The  second  term  for  the  expressions  of  and  G£  in  Equation 

21  is  the  residue  contribution  of  the  poles.  For  an  earth  model  with 

no  anelastic  attenuation  loss,  the  poles  lie  on  the  real  axis  in  the 

range  k < k < ka  The  poles  are  found  by  locating  the  zeros  of 
eN  - Br 

the  function  F.  by  a search  in  the  region  k.  + Ak  < k < k.  to  find 

L BN  ~ “ 61 

zero  crossings  of  tne  function,  after  which  an  interval  halving  process 

is  used  to  refine  the  root.  To  avoid  numerical  problems  at  the  branch 

point,  a possible  pole  at  k * k0  is  neglected.  This  is  a minor  error 

if  Ak  is  small  enough.  A value  Ak  = (k.  - kQ  )/100  is  used  in  the 

B1  bN 

program.  For  reasonable  continental  crustal  models  with  crustal  thick- 
nesses of  about  40  km,  this  choice  of  Ak  should  be  adequate  to  frequen- 
cies of  3 Hz,  since  the  number  of  poles  increases  with  increasing  fre- 
quency. The  residue  contribution  of  a function  g(z)/h(z),  having 
a zero  of  h(z)  at  z * zQ,  is  very  simply 

9<V/h',Z0>* 

where  h'(z)  is  the  first  derivative  of  h(z). 

30.  The  third  term  in  the  expressions  for  G^  and  G2  is  an  integral 
of  the  form 

/ fO’x.wjK^xrJx  dx,  (25) 

where  f(ix,u>)  is  a complex  quantity.  The  function  Kn(z)  decreases  in 
an  exponential  manner  for  Increasing  z and  has  a singularity  at  z * 0. 

The  exponential  decrease  with  large  values  of  the  argument  suggests  the 


1 


i 

use  of  a Gauss -Laguerre  Integration  rule  If  f(i-r,w)  does  not  oscillate 
too  fast.  After  a change  of  variable  and  application  of  the  Gauss- 
Laguerre  Integration  rule.  Equation  25  Is  approximated  by 

i m jy 

~2  ^ wk  [ f(^k,u)  xk  exp{xk)Kn(xk\ )]  (26) 

where  xk  and  wk  are  the  abscissas  and  weights  of  an  order  Gauss- 
Laguerre  rule. 

31.  Error  is  Introduced  Into  the  evaluation  of  Equation  26 

because  of  the  oscillatory  nature  of  f(ix,u).  This  can  be  mitigated 

by  using  a very  high  order  rule,  so  that  the  abscissas  are  closely 

enough  spaced.  Because  the  weights  wk  decrease  very  fast  for  large 

values  of  the  Index  k,  one  does  not  need  to  use  all  m terms  Of  the 

summation.  The  computer  program  WESHASK  uses  the  first  24  weights 

and  aoscissas  of  a m = 68  order  rule  which  were  taken  from  the 

9 2 

tables  of  Stroud  and  Secrest.  Because  of  the  r term  in  the 
expression,  the  contribution  of  Equation  26  to  the  total  solution 
only  becomes  Important  at  short  distances.  The  approximation  used 
in  the  computer  program  gives  good  results  for  distances  greater  than 
5 km  and  for  frequencies  less  than  1 Hz. 

Independent  Solution  for  Simplified  Case 

32.  To  test  the  computer  program  WESHASK,  an  independent 
solution  was  found  for  the  simple  case  of  a source  within  a single 
layer  overlying  a halfspace.  For  a source  at  a depth  h within  a 
layer  with  thickness  d^,  the  expressions  for  G^  and  Gg  of  Equation  11 
take  the  form 


k JQ(kr)  dk  (27a) 


G,  = / 

il  O 


3?  • P2  , 

v S + — v C 
-T  3,  P-,  v02  B, 

J2 

8] 

p*  V“1 + 
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“here  se,  a s6)<dl  * hl>-  S,  * S,(dt  ' V-  % ■ % <d, ) 


and 


C3^  = Cg  otlier  symbols  are  as  previously  defined 

33.  Defining  the  1 
reflection  coefficient  as 


33.  Defining  the  layer  rigidity  as  nn  * Pnen  and  a complex 


- R^exp[-v^ (4dj-h)J  + R2exp[-vg^ (4d]+h)]  +...+ 


(28) 


-k 


0 plV 


exp[-v  hj  - R exp[-v  (2d,  - 
■}  el  1 


h)3 


- R exp[-vg^(2d^+h)]  + R2exp[-v^(4d1-h)] 


+ R2exp[-v0  (4d,+h)]  +...+ 
P1  1 


lb 
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G2  = 


ef  . p2%  , 


f 

0 


p2Vb,  * C17vB,S6 

*2 


e 


k^  Jj(kr)  dk 


1 


2 P1  P1  P1 


(27b) 


p2V  ‘ 

R * £ 1 , 

V2\  + Va, 

Equation  27  can  be  rewritten  as  an  infinite  series  as  follows: 


«»  i 

G,  = / — 

• 0 P] 


exp(-v^h)  + Rexp[-vg^  (2d]-h)j  -R  exp[-v^  (2d]+h)] 


34.  This  expansion  is  termed  a generalized  ray  expansion  since 
each  term  represents  the  contribution  of  rays  that  are  multiply  re- 
flected after  leaving  the  source.  The  first  term  in  the  Gj  and 
expansions  represent  the  contribution  of  the  direct  ray  from  the 
source  to  the  receiver  at  the  free  surface,  while  the  next  two  terms 
represent  rays  that  are  reflected  or  refracted  once  from  the  inter- 

O 

face  at  a depth  z * dj.  The  term  generalized  ray  is  used  since  the 
integral  of  each  term  yields  both  the  reflected  and  refracted  arrivals 
associated  with  each  ray.  In  an  infinite  medium,  just  the  first  term 
comprises  the  solution  for  G^  and  G2  since  R * 0. 

35.  The  first  terms  of  the  expressions  for  Gj  and  G2  in 
Equation  28  can  be  evaluated  directly  by  taking  derivatives  of  the 
Sommerfeld  integral.  It  can  be  easily  shown  that 

/ exp[-vgh]Jo(kr)k  dk  = £ [ exp(-1kgR)  = H] 

/ ^ exp[-v  h]J.(kr)k  dk  = ^ f ^ exp(-ik  R)  = H- 

0 6 1 R t $R  R*J  6 2 

2 2 2 

where  R = h+  r. 

36.  Equation  29  can  be  used  to  evaluate  the  individual  terms 
of  Equation  28.  The  two  basic  integrals  to  be  evaluated  are  of  the 
form 

00 

/ exp[-vfl  Z]  Jn(kr)k  dk  (30a) 

o p-j  u 

and 

/ ~ r"  exp[-vfi  Z]  J,(kr)  k dk  . (30b) 

0 V81  B1  1 

The  evaluation  of  these  two  integrals  is  complicated  by  the  oscilla- 
tory character  of  the  Bessel  function  for  large  k.  However,  for 
large  k,  the  complex  reflection  coefficient  R approaches  a limit. 

Mo  - Mi 

R(“.«)  ■ r—  • (31 ) 

m2  + 


(29) 

I 
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37.  The  oscillatory  behavior  for  large  k can  have  a minimal  effect 
if  the  integrals  of  Equation  30  are  rewritten  using  Equations  29  and  31 
in  the  following  manner: 

/ [ ^(k.uO  - rV.u)]  exp[-v  Z]  Jn(kr)k  dk  + rV.uOH,  (32a) 

0 P-j  u I 

/ J_  [ Rn ( k ,0) ) - R?(-,«)3  expC-v.  Z]  J, (kr)k  dk  + Rn(»,u))H?  . (32b) 

o vg^  61  I ^ 

38.  The  integrals  to  be  evaluated  in  Equation  32  are  of  the 


form 

and 


oo 

Ii  s / g(k)exp[-v  Z]J  (kr)k  dk 

• o p-j  o 

U * / — - g(k)exp[-v  Z]J,(kr)k  dk, 
i 0 B-|  61  1 


(33a) 

(33b) 


where  g(k)  is  function  which  goes  to  zero  for  large  k. 

39.  A critical  point  in  the  evaluation  of  Equation  33  is  at 

k = k where  vQ  changes  from  positive  imaginary  to  positive  real 
el  el 

with  increasing  k.  To  properly  take  into  account  the  change  in 
character  of  the  integrands,  the  integration  is  performed  over 
two  ranges  oy  using  the  following  transformations: 
a.  0 < k < k 


k * k_  siny 
B1 

dk  * k.  cosy  dy 
61 


V.  = ik.  cosy  , for  0 < Y < n/2 
61  B1  - - 

b.  kc  < k < « 

- B1  - 

k * k coshn 
61 

dk  = k.  sinhn  dn 
61 
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v-  = k0  sinh  n,  for  0 < n < «°  . 
el  61 

40.  Using  these  transformations.  Equations  33a  and  33b  become 
the  following: 


tt/2 


I,  = / g(k  sinY)exp[-ik  ZcosY]Jn(k  r sinY)k  cosYsinY  dY 

I 0 P]  u »i 


2 

+ / g(k„  coshn)exp[-k.  Zs1nhn]Jn(k0  rcoshn)k0  coshnsinhndn 

OP1]  p^j  u ti'j  ^1 

(34a) 


tt/2 


I,  » -i  / g(k«  sinY)exp[-ik  ZcosY]J,(k  r sinY)k  sin  Y dY 
c.  0 p]  P]  ^1  Pi 


2 2 

+ / g(kfl  coshn)exp[-ka  Zsinhn]vMk0  rcoshn)k„  cosh  n dn 

0 Pi  P]  * ^1  ^1 

(341 

41.  The  first  integral  in  the  expressions  for  1^  and  is 
easily  evaluated  using  a trapezoidal  rule.  The  limitations  on  the 
applicability  of  the  trapezoidal  rule  are  as  discussed  previously  in 
Paragraph  28.  The  second  integral  in  each  expression  can  be  evalu- 
ated using  the  Gauss-Laguerre  integration  rule  after  one  last  trans- 
formation. Let  x = k Z sinh  n.  Then 

pi 


and 


dx  = kQ  Z cosh  n dn 

61 

[k«  + (x/Z)2]1/2  = k cosh  n. 

61  61 

42.  The  numerical  approximations  to  ^ and  I2  are  finally 

M-l  ? 

i!  - I g[k(Ym)]  exp[-ik  Zcos  Ym]J  [k(Ym)r]k‘  sin  Ym  cos  ym  ay 
m =1  1 1 

+ \ I wj  g[k(Xj)]  Xj  JQ[k(Xj)r]  (35a) 

Z j"  • 
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(35b) 


‘2  ' _,il  SCk(Tn,):i  exp[-ilt61Zcos  ’WI¥lt,’f.i,r:|lt«,  Sln2vm  4y 

+ - l W.  g[k(x.)]  k(x-)  JJMxJr]  , 
z j=i  j j j i a 

where 

Ay  = ir/2M 

rm  • » 4y 
HyJ  = kB)sin  Ym 

and  Xj  and  are  the  abscissas  and  weights  of  an  n-point  Gauss- 
Laguerre  integration  rule  and 

k(x.)  = [k2^  + (x./Z)2  ]1/2  * 

Tne  techniques  outlined  in  this  section  are  used  in  the  computer 
program  WESREFt  which  is  given  in  Appendix  B. 

Numerical  Tests 

43.  In  the  process  of  generating  realiatic  seismograms,  a 
fast  Fourier  transform],  FFT,  is  used  to  perform  the  integration  re- 
quired by  Equation  12  to  convert  the  frequency  domain  representation 
of  the  solution  into  the  desired  time  series.  One  problem  that  arises 
is  the  representation  of  a unit  step  function  by  the  FFT.  This  arises 
because  the  FFT,  unlike  the  Fourier  transform,  is  periodic  in  time.^ 

To  avoid  numerical  problems  associated  with  discontinuities,  the 
velocity  on  the  fault  is  specified  rather  than  the  displacement.  This 
means  that  the  seismograms  computed  are  ground  motion  velocity  histories 
rather  than  ground  motion  displacement  histories.  The  solution,  though, 
is  amenable  to  either  integration  or  differentiation  to  form  displace- 
ment or  acceleration  time  histories. 

44.  The  source  function  used  to  represent  the  velocity  of  the 
motion  of  the  fault  is  given  by  the  function  s(t),  defined  as 
follows: 


2xs(t) 


0 

t < 0 

(t/x)2 

0 < t < x 

(t/t)2  + 2 ( t/x ) - 1 

x < t < 3x  (36) 

\ (t/x)2  - 4 ( t/x ) + 8 

3x  < t < 4x 

0 

t > 4x 

The  pulse  s(t)  is  defined  such  that  the  area  under  the  pulse  is 
equal  to  one. 

45.  For  the  examples  that  follow  the  values  t = 0.5,  1,  2 sec 
have  Deen  used.  To  avoid  numerical  problems  at  high  frequencies, 
s(t)  is  low  pass  filtered  to  pass  frequencies  less  than  1.0  Hz. 
Figures  2 and  3 show  the  low  pass  filtered  pulse  and  its  Fourier 
amplitude  spectrum  for  x = 0.5  sec  and  x “ 1.0  sec,  respectively. 

The  cutoff  frequency  of  1.0  Hz  was  chosen  to  coincide  with  a spectral 
minimum  of  S(f)  for  x = 0.5,  1 , or  2 sec.  This  choice  prevents  sharp 
discontinuities  in  the  frequency  domain,  which  would  introduce  noise 
in  the  time  domain  signal.  This  technique  eliminates  some  of  the 
high  frequency  noise  problem  that  was  present  in  earlier  work 

by  Herrmann  and  NuttliJ^2 

46.  The  earth  model  used  to  test  the  computer  program  WESHASK 
by  comparison  to  the  results  of  the  program  WESREFL  is  a simplified 
continental  model  (SCM)  given  in  Table  1.  This  model  represents  a 
first  approximation  to  a realistic  earth  model  for  the  central 
United  States.  A second,  more  realistic  earth  model  is  also  given 
in  Table  1 under  the  heading  Central  United  States  Model  (CUS). 

This  second  model  is  in  reasonable  agreement  with  P wave  refraction 
and  surface  wave  dispersion  models  proposed  for  the  central  United 
States.'3-'4-'5-'6 

47.  As  input  to  the  computer  programs  WESHASK  and  WESREFL, 


25 


RMP  (CM-SEC) 


in 

o 


J\ 

-fe.on  ib. do  2b.  oo  «b.oo  sb.oo  ii. 


FREQ  (HZ) 

Fig.  3.  Filtered  source  pulse  and  amplitude  spectrum  for  t = 1.0  sec. 
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Layer 


Thickness 

P Vel 

S Vel 

Density 

d km 

a km/sec 

B km/sec 

p gm/cm* 

Simplified  Continental  Model  (SCM): 


40 

6.15 

3.55 

2.8 

— 

8.09 

4.67 

3.3 

Central  United  States  Model  (CUS) 

1 

5.00 

2.89 

2.5 

9 

6.10 

3.52 

2.7 

10 

6.40 

3.70 

2.9 

20 

6.70 

3.87 

3.0 

-- 

8.15 

4.70 

3.4 

► 


layer  thicknesses  and  distances  are  given  in  units  of  km,  the  layer  wave 
velocities  are  given  in  units  of  km/sec  and  the  density  is  given  in 

3 

units  of  gm/cm  . These  units  were  chosen  for  ease  of  input  of  numerical 
data.  In  line  with  this  usage,  the  ground  motion  velocity  time  histo- 
ries are  in  units  of  cm/sec  if  it  is  understood  that  a source  of  seismic 
20 

moment  10  dyne-cm  is  used.  To  compute  the  ground  motion  for  a seismic 
24 

moment  of  10  dyne-cm,  for  example,  one  must  multiply  the  computer  out- 

4 

put  by  the  factor  10  . Note  also  that  the  output  given  is  the  function 
Gj  or  G2  and  the  results  must  be  multiplied  by  the  functions  F^  or  F2  of 
Equation  11  to  give  the  proper  results  for  comparison  to  real  data. 

48.  Figure  4 is  the  comparison  test  of  the  G^  solution  at  a dis- 
tance r = 10  km  from  a source  10  km  deep  in  the  SCM  model.  ’A  source 
with  t = 1 sec  is  used.  To  understand  the  relative  effect  and  impor- 
tance of  each  term  in  Equation  21,  the  contribution  of  each  additional 
term  to  the  final  solution  is  shown.  For  ease  of  presentation,  the 
time  series  is  plotted  as  a function  of  the  reduced  travel  time 
"t  - x/4.57"  where  t is  the  true  travel  time  and  x is  the  distance 
from  the  source.  Figure  4a  shows  the  pole  contribution  to  the  total 
solution.  Figure  4b  shows  the  effect  of  adding  the  contribution  of 
the  branch  line  integral  along  the  real  k axis  to  the  pole  contribu- 
tion. Figure  4c  shows  the  complete  solution  provided  by  WESHASK, 
e.g.,  that  which  includes  the  pole  contribution  and  the  contributions 
of  the  branch  line  integrals  along  the  real  k axis  and  along  the  neg- 
ative imaginary  k axis.  Figure  4d  is  the  solution  provided  by  WESREFL. 
It  is  seen  that  the  pole  contribution  is  non-causal  and  of  low  ampli- 
tude. (The  later  arrivals  are  in  fact  early  negative  time  arrivals 
because  of  the  periodicity  of  the  FFT).  The  addition  of  the  branch 
line  integral  along  the  real  axis  improves  causality  and  raises  the 
signal  amplitude  to  the  final  level.  The  branch  line  integral  along 
the  negative  k axis  affects  the  amplitude  level  only  slightly  while 
making  the  resultant  signal  causal.  The  agreement  between  the  WESHASK 
(Figure  4c)  and  the  WESREFL  (Figure  4d)  solutions  is  excellent.  The 
slight  motion  at  a reduced  travel  time  of  about  15-20  sec  is  due  to 
energy  reflected  once  from  the  layer  boundary  at  a depth  of  40  km. 
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a.  G1  solution  at  10  km,  pole  contribution  only; 

b.  Pole  contribution  and  real  axis  branch  line  Integral; 

c.  Pole  contribution  and  both  branch  line  Integrals; 

d.  Ray  theory  solution 
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Figure  5 is  similar  to  Figure  4 except  that  the  G2  solution  is  shown  for 
the  same  source  and  receiver  condit’‘'ns. 

49.  Figure  6 compares  the  G-j  solution  obtained  at  a distance  of 
300  km  from  the  same  source  used  to  generate  Figure  4.  This  figure 
shows  that  at  large  distances  the  pole  contribution  (Figure  6a)  is  the 
most  important  contributor  to  the  final  solution.  The  addition  of  the 
branch  line  integral  along  the  real  axis  to  the  pole  contribution  (Fig- 
ure 6b)  makes  the  signal  causal  (e.g.,  no  arrivals  before  predicted  Ray 
theory  arrivals).  As  expected  from  Equation  26,  the  branch  line  inte- 
gral along  the  negative  imaginary  k axis  makes  an  insignificant  contri- 

_2 

bution  at  large  distances  because  of  its  r character.  Again  the  com- 
parison between  the  final  solution  obtained  using  WESHASK  (Figure  6c) 
and  the  WESREFL  solution  (Figure  6d)  is  very  good.  Figure  7 is  the  G2 
solution  corresponding  to  the  G.|  solution  of  Figure  6. 

50.  While  no  independent  has  been  developed  yet  to  test  WESHASK 
for  more  complicated  earth  models  than  a single  layer  overlying  a half- 
space, it  is  of  value  to  consider  the  importance  of  the  various  contri- 
butions of  the  terms  of  Equation  21.  Figure  8 is  the  G^  solution  for 

a source  at  a depth  of  10  km  in  the  CUS  earth  model.  A source  pulse 
with  t = 0.5  sec  is  used.  From  top  to  bottom,  the  various  traces  are 
the  pole  contribution,  the  contribution  of  the  poles  and  real  axis 
branch  line  integral,  and  the  complete  solution.  As  for  the  simple 
earth  model  at  short  distances,  the  pole  contribution  alone  is  non- 
causal  and  of  low  amplitude.  The  addition  of  the  branch  line  integral 
along  the  real  k axis  improves  causality  and  raises  the  amplitude  to 
its  final  level.  The  addition  of  the  branch  line  integral  along  the 
negative  imaginary  k axis  makes  the  signal  causal.  Figure  9 is  the  G2 
counterpart  of  Figure  8.  The  slight  overshoot  in  the  main  pulse  is  an 
effect  of  only  using  frequencies  less  than  1.0  Hz. 

51.  Figure  10  is  the  G^  solution  at  a distance  of  10  km  from  a 
source  at  a depth  of  10  km  in  the  CUS  model.  In  this  case  a longer 
source  pulse  was  used  with  t * 1 sec.  Figure  11  is  the  G2  solution 
corresponding  to  the  G^  solution  of  Figure  10.  A comparison  of 
Figures  8-11  shows  that  the  causality  and  amplitude  errors 
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Fig.  5.  a.  G.,  solution  at  10  km,  pole  contribution  only; 

b^.  Pole  contribution  and  real  axis  branch  line  integral: 
c^  Pole  contribution  and  both  branch  line  Integrals; 
d.  Ray  theory  solution 
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Fig.  6.  a.  G,  solution  at  300  km,  pole  contribution  only; 

■ 

d.  Pole  contribution  and  real  axis  branch  line  integral; 
c^  Pole  contribution  and  both  branch  line  Integrals; 
d.  Ray  theory  solution 


I a-  !<■> 


T - X/4.67 


Fig.  7. 


a.  G2  solution  at  300  km,  pole  contribution  only; 

b.  Pole  contribution  and  real  axis  branch  line  Integral; 
, Pole  contribution  and  both  branch  line  Integrals; 

, Ray  theory  solution 
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Fig.  8.  Comparison  of  contribution  to  Gj  for  CUS  model  at  10  km.- 
From  top  to  bottom:  pole  contribution;  pole  contribution 
plus  real  axis  branch  line  integral;  complete  solution. 
t = 0.5  sec. 
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Fig.  9.  Comparison  of  contribution  to  G2  for  CUS  model  at  10  km. 

From  top  to  bottom:  pole  contriDution;  pole  contribution 
plus  real  axis  branch  line  integral:  complete  solution. 
t «=  0.5  sec. 


Fig.  10.  Comparison  of  contribution  to  for  CUS  model  at  10  km. 
From  top  to  bottom:  pole  contribution;  pole  contribution 
plus  real  axis  branch  line  integral;  complete  solution, 
x * 1 .0  sec*. 
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Fig.  11.  Comparison  of  contribution  to  for  CUS  model  at  10  km. 
From  top  to  bottom:  pole  contribution;  pole  contribution 
plus  real  axis  branch  line  integral;  complete  solution, 
r = 1 .0  sec. 
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associated  with  using  a partial  solution  obtained  by  ignoring 
one  term  of  the  complete  expansion  of  Equation  21  becomes  worse 
as  a longer  source  pulse  Is  used.  Figures  12  and  13  are  the  Gj  and 
comparisons  at  a distance  of  50  km  from  a source  at  a depth  of 
10  km  in  the  CUS  model  with  a pulse  having  t * 1 sec.  The  pole  con- 
tribution by  itself  yields  amplitudes  of  the  right  order.  The  Gj 
solution  requires  the  branch  line  integrals  for  causality,  while  the 
Gj  solution  is  almost  completely  described  by  the  pole  contribution 
al one . 

52.  The  value  of  the  exercise  just  performed  by  studying  the 
solutions  shown  in  Figures  8 - 13,  is  that  distance  ranges  at  which 
computational  efficiency  can  be  obtained  are  indicated.  For  example, 
at  short  distances,  the  complete  solution  complete  with  poles  and 
branch  line  integrals  must  be  obtained.  At  large  distances,  especially, 
if  one  is  not  interested  in  causality,  the  pole  contribution  by  Itself 
suffices  to  provide  a very  realistic  estimate  of  the  final  solution. 

By  not  having  to  perform  the  branch  line  integral  along  the  negative 
imaginary  axis  at  distances  greater  than  100  km,  for  example,  con- 
siderable computer  time  can  be  saved  in  obtaining  the  final  solution. 

53.  Figure  14  shows  the  dependence  of  G^ , left  side,  and  G2, 
right  side,  upon  distance.  The  source  is  at  a depth  of  10  km  in 
the  SCM  model.  A source  pulse  with  t * 1 sec  is  used.  Secondary 
arrivals  due  to  deep  reflections  are  very  prominent  in  the  G^  solution 
since  the  vertical  dip  slip  source  is  very  efficient  at  generating 
waves  that  leave  the  source  In  a near  vertical  direction.  On  the 
other  hand  the  vertical  strike  slip  source  or  45°  dip  slip  source 

is  not  very  efficient  at  generating  near  vertical  waves.  Hence, 
the  near  vertical  reflections  are  more  prominent  in  the  Gj  solution. 

At  larger  distances,  rays  which  left  the  source  near  vertically  at 
short  distances  reflect  off  the  interface  at  the  depth  of  40  km  with 
angles  near  the  critical  angle.  As  a ray  becomes  supercritical ly 
reflected,  there  is  little  energy  loss  upon  refiction  and  a phase 
change  Is  Introduced.  It  is  seen  from  Figure  14  that  very  low  ampli- 
tude phases  at  short  distances  become  large  contributors  to  the  signal 


Fig.  14.  Comparison  of  G-j  solution,  a,  and  G2  solution,  b, 
for  SCH  model  as  a function  of  distance  In  km 


from  a source  at  depth  of  10  km  with  x s 1.0  sec_ 
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at  large  distances  for  precisely  this  reason.  It  is  also  seen  how 
the  various  multiple  reflections  combine  together  to  form  what  is 
known  as  a surface  wave  at  large  distances.  The  refraction  arrival 
is  of  relatively  small  amplitude  compared  to  the  surface  wave  at 
large  distances.  The  progression  of  the  SH  wave  signal  from  a simple 
pulse  at  short  distances  to  the  surface  wave,  called  Love  wave,  is 
gradual.  The  surface  wave  can  be  said  to  take  over  from  the  direct 
pulse  at  a distance  at  which  the  first  supercritically  reflected 
waves  arrive.  At  tnis  point  the  waveform  changes  from  a simple  pulse 
to  a complicated  set  of  arrivals. 

54.  Figure  15  compares  the  6^  and  G2  solutions  as  a function 
of  distance  for  a source  with  t = 1.0  sec  at  a depth  of  10  km  in  the 
CUS  earth  model.  There  seems  to  very  little  difference  between  this 
solution  and  that  presented  in  Figure  14  at  short  distances.  However, 
at  large  distances,  the  arrivals  are  no  longer  as  distinct  because 
trie  more  complicated  earth  model  permits  many  more  reflections  and 
refractions  by  the  various  layers.  Figure  16  shows  the  G1  and  G2 
solutions  as  a function  of  distance  for  a source  with  t * 0.5  sec 
at  a depth  of  10  km  in  the  CUS  earth  model.  With  the  shorter  source 
pulse  the  various  crustal  reflections  can  be  easily  seen  compared 
to  Figure  15  with  tne  longer  source  pulse.  A close  examination  of 
the  records  indicates  that  numerical  noise  is  present  only  at  the 
very  short  and  very  large  distances,  as  is  expected  from  the  discussion 
of  numerical  Integration  techniques  used. 


Fig.  15.  Comparison  of  G^  solution,  a^  and  Gg  solution,  b^, 
for  CUS  model  as  a function  of  distance  In  km 


from  a source  at  a depth  of  10  km  with  x ■ 1.0  sec 


PART  IV:  SUMMARY  AND  CONCLUSIONS 


55.  The  numerical  results  presented  in- the  previous  section 
show  that  the  tasks  originally  undertaken  for  this  research  project 
have  been  accomplished.  For  the  first  time  a workable  computer  pro- 
gram package  is  available  for  making  realistic  predictions  of  SH  wave 
ground  motion  both  near  the  seismic  source  and  at  large  distances 
from  the  source.  The  solutions  presented  seem  valid  in  the  distance 
range  of  5 - 500  km  for  frequencies  less  than  1 Hz.  Future  work 
should  be  directed  toward  the  extension  of  the  computer  program  to 
shorter  epi central  distances  and  higher  frequencies. 

56.  This  ground  motion  model  should  be  of  use  in  its  present 
form  to  make  some  deductions  as  to  the  form  of  an  empirical  model  for 
fitting  observed  ground  motion.  Many  attempts  at  modeling  strong 
motion  data  have  been  made  by  trying  to  fit  observed  time  histories 

to  motions  expected  from  a dislocation  source  in  an  infinite  medium.^ 
This  technique  ignores  the  generation  of  surface  waves  and  the  complex- 
ity introduced  by  the  free  surface.  It  is  equivalent  to  using  only 
the  first  term  of  the  6^  and  G,,  expansions  of  Equation  28.  The  adequacy 
of  this  elementary  modeling  can  be  tested  using  the  computer  programs 
developed  in  this  report. 

57.  Figure  17  demonstrates  the  effect  of  variation  of  focal 
depth  upon  the  maximum  ground  velocity  for  a source  with  t * 1 sec 
at  depths  of  5,  10  and  20  km  in  the  SCM  earth  model  of  Table  1.  The 
velocity  values  must  be  adjusted  as  described  in  Paragraph  47  before 
comparing  to  real  data.  The  solid  curves  indicate  the  theoretical 
solutions  for  the  halfspace  having  the  material  properties  of  the 
first  layer  of  the  SCM  earth  model.  It  Is  seen  that  the  halfspace 
solution  is  adequate  out  to  about  75  km.  This  is  because  at  75  km 
for  this  earth  model  the  first  supercritical  reflections  occur  and 
the  surface  wave  begins  to  form.  The  difference  in  the  geometrical 
spreading  between  the  Gj  and  solutions  is  due  to  the  difference 
in  radiation  patterns  between  the  two  basic  solutions.  Figure  18 
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17.  Effect  of  focal  depth  variation  (5,  10,  20  km)  udoo  maximum  ground 
as  a function  of  distance  for  t * 1.0  sec 
using  SCM  model 


Fig.  18.  Effect  of  source  pulse  length  (t  = 0.5,  1,  z sec;  upon  maximum  ground 
velocity  as  a function  of  distance  for  depth  of  10  km  in  SCM  model 


shows  the  effect  of  variation  of  source  pulse  length  for  a source  at 
a depth  of  10  km  in  the  SCM  earth  model.  The  maximum  ground  velocity 
was  computed  using  t = 0.5,  1 and  2 sec.  For  a constant  area  under 
the  source  pulse,  or  in  other  terms  a constant  seismic  moment,  the 
maximum  ground  velocity  increases  as  t decreases.  The  effect  of 
the  transition  from  a pulselike  arrival  to  a surface  wave  is  again 
seen  at  a distance  of  75  km  by  the  departure  of  the  computed  solution 
for  the  SCM  earth  model  indicated  by  the  symbols  and  the  halfspace 
approximation  given  by  the  solid  lines. 

58.  Figures  17  and  18  were  obtained  for  a very  simple  earth 
model.  The  computer  program  WESHASK  can  be  used  to  do  a similar 
study  for  other  earth  models  and  source  parameters.  For  a crustal 
model  with  a shallower  depth  to  the  crust-mantle  interface,  simple 
model  scaling  indicates  that  the  transition  from  a near-field 
pulselike  ground  motion  character  to  a far-field  surface  wave 


character  would  occur  at  shorter  epicentral  distances  than  for  the 
SCH  earth  model.  It  is  quite  feasible  to  determine  figures  similar 
to  Figures  17  and  18  for  a southern  California  earth  model  in  order 


to  obtain  an  insight  as  to  the  empirical  shape  to  use  when  fitting 
southern  California  strong  motion  data  as  a function  of  distance. 

59.  Aside  from  giving  an  insight  into  the  nature 
of  gross  ground  motion  parameter  modeling,  such  as  maximum 
velocity  versus  distance,  the  computer  program  developed  here 
should  be  able  to  predict  strong  motion  ground  displacement 
time  histories  for  real  earthquakes  as  well  as  the  generalized 
ray  technique  of  Helmberger  and  Malone.1  To  predict  strong 
motion  acceleration  time  histories,  whose  high  frequency 
content  is  strongly  affected  by  the  heterogeneities  of  the 
real  earth,  further  research  is  required  on  the  statistical 
modeling  of  high  frequency  waves  scattered  from  these 
heterogeneities.  When  this  is  done,  a technique  proposed 
by  Herrmann18  could  be  used  to  yield  a synthetic  seismogram 
having  the  coherent  low  frequency  character  described  by 
the  theory  presented  in  this  report  and  the  high  frequency 
content  of  the  scattered  waves.  Such  a ground  motion  time 
history  would  be  as  realistic  as  is  presently  possible. 
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